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Solvation of small and large clusters are studied by simulation, considering a range of solvent-solute 
attractive energy strengths. Over a wide range of conditions, both for solvation in the Lennard- Jones 
liquid and in the SPC model of water, it is shown that the mean solvent density varies linearly with 
changes in solvent-solute adhesion or attractive energy strength. This behavior is understood from 
the perspective of Weeks' theory of solvation [Ann. Rev. Phys. Chem. 2002, 53, 533] and supports 
theories based upon that perspective. 



I. INTRODUCTION 

Some solute species dissolve in a liquid solvent over 
a large range of conditions. Others demix or phase 
separate from the liquid solvent even at relatively low 
solute concentrations. This paper focuses on meso- 
scopic and microscopic manifestations of these macro- 
scopic phenomena. Our motivation is to understand 
the range of solvation behaviors that are possible in 
aqueous solutions when hydrophobic species are al- 
tered to hydrophilic and vice versa. While the mo- 
tivation concerns this specific class of solutions, the 
physical principles involved are relatively general. 

The basic theory has been formulated by Weeks [1] . 
One first considers the effect of the solute in exclud- 
ing a region of space from the solute. The average 
effect is described by the response of the solvent to 
a so-called unbalancing potential. Unbalancing arises 
because the excluded volume presents a region where 
solvent molecules lose adhesive forces. In a homoge- 
neous fluid, these adhesive forces tend to be in bal- 
ance, often allowing their effects to be described as a 
uniform field or chemical potential shift [2J |3] . Intro- 
ducing a large enough excluded volume, however, de- 
stroys this homogeneity, and the response can be very 
large, particularly so if the liquid exists near coexis- 
tence with its vapor. Indeed, in the macroscopic limit, 
the average liquid density adjacent to excluded regions 
bounded by a plane is /3p, where p is the pressure, 
and I3~^ = k^T is temperature times Boltzmann's 
constant [1]. In this limit, therefore, the mean liquid 
density at the surface is a bulk property. It exhibits a 
phase transition reflecting the bulk liquid- vapor phase 
transition. Weeks' unbalancing potential captures the 
effects of this transition in a mean molecular field ap- 
proximation. 

A typical solute does not exclude liquid from a 
macroscopic volume. Nevertheless, the unbalancing 
potential can be substantial and the response to it 
significant, though not singular, provided the solvent- 
exposed surface area of the solute has low enough cur- 
vature and large enough size. Weeks' expression for 



the unbalancing potential shows that its strength is 
roughly proportional to the surface area, A. Thus, to 
linear order in the response, 

{p{v)),^a{v)-Ah{v). (1) 

Here, p{t) denotes the solvent density field at position 
r, which we take to be outside the excluded volume, 
and the pointed brackets indicate equilibrium ensem- 
ble average, with the subscript implying that only 
excluded volume forces act between solute and solvent. 
The function a(r) and 6(r) are properties of the sol- 
vent and the solute position and shape; a(r) gives the 
response in the absence of unbalancing forces; A b{r) is 
significant in size only when the liquid solvent is close 
to phase coexistence with its vapor and A is compa- 
rable to or larger than the surface area of the critical 
nucleus for the phase transition; for positions r that 
are beyond a bulk correlation length of the solvent, 
a(r) approaches the bulk liquid density, and b(r) ap- 
proaches zero. 

Fluctuations around the average and the effects of 
forces beyond those of excluded volume are captured 
in the second step of Weeks' treatment. In this step, 
one notes that solvent density fluctuations in a ho- 
mogeneous hquid obey Gaussian statistics [S] [HI [7], 
and assumes this statistics is more generally valid for 
all small length scale fluctuations in a liquid, even in 
the presence of inhomogeneities. Therefore, if the so- 
lute attracts the solvent molecules with a potential 
strength proportional to A, the effect of the attrac- 
tion on the mean solvent density is linear in A, i.e., 
(p(r))^ » (p(r))o + •^'S'(r) , where S{r) is the net or 
integrated response of the solvent density fleld at r to 
the attractions from the solute. If these attractions 
emanate from interaction sites in the solute and are 
flnite ranged, then S(r) will be approximately linear 
in A, leading to 

{p{v)),»{p{r)), + XAs{r), (2) 

where s(r) is a property of the solvent and the func- 
tional form of the attractive interaction. Finally, as 
both 6(r) and s(r) describe responses of the solvent 
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near the solute and per unit exposed area of the solute, 
one may expect that integrated or coarse-grained ef- 
fects of the two are approximately proportional. That 
is to say, 6(r) w A*s(r), where the over-bar indicates 
an integral or coarse-graining over some specified mi- 
croscopic length scale. The value of the proportion- 
ality constant. A*, depends upon the specific solvent, 
the specific solute-solvent interactions, and the spe- 
cific coarse-graining perscription. Adopting this pro- 
portionality gives 



(p(r))^«a(r) + (A-A*) As{v). 



(3) 



Equations ([iJ, ^ and ^ are not generally true, as 
they do not properly account for saturation in either 
the limit of very large A or the limit of very large A. 
The point of this paper, however, is that these equa- 
tions are good approximations in a variety of physi- 
cally pertinent circumstances. We demonstrate their 
accuracy in the next section by presenting results of 
computer simulation calculations of solvent density 
fields and solvation free energies for solute clusters 
of various sizes and solvent-solute attractive energy 
strengths, considering both the Lennard- Jones liquid 
solvent and water solvent. The results indicate that 
variations from favorable to unfavorable solvation, i.e., 
the changes in solvation in passing from large A to 
small A, occur smoothly and in ways that are easy to 
anticipate. Further, the results provide confidence in 
theories built from Weeks' perspective. These impli- 
cations are discussed in Section III. 



II. SIMULATION RESULTS 

To elucidate the response of a liquid to the presence 
of a solute of general shape and interaction strength, 
we consider a series of models for both solutes and 
solvent of increasing complexity. 



A. Spherical solutes in Lennard— Jones solvent 

First we consider the case of a spherical solute im- 
mersed in the Lennard- Jones solvent, i.e., a fiuid of 
particles that interact through the pair potential 
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(4) 



In this equation r denotes the distance between two 
solvent particles, and e and a define the energy and 
length scales of the simulation. The solvent-solvent 
potentials are shifted so as to be continuous and zero 
beyond the cut-off distance Tc = 2.5cr. The reduced 
thermodynamic state conditions for the solvent are 



fegT/e = 0.85, pa^/e = 0.022 and pa^ = 0.7, where 
p is the bulk density. These conditions place this 
Lennard-Jones liquid solvent at or very close to liquid- 
vapor coexistence [8J. 

The solute is a sphere of radius R, and its inter- 
action with the solvent is that of a sphere consisting 
of a Lennard-Jones material with homogenous den- 
sity equal to that of the solvent. A solvent particle 
at distance r from the solute center has the potential 
energy 13 
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where r± = r ± R, ii r > R, and (j){r) = oo otherwise. 
This potential is minimal at a distance r-aan, and its 
minimum value <i>niin = *i'(^min) depends on the solute 
size R. 

To systematically investigate the role of attractive 
interactions, we employ the scheme introduced by 
Weeks, Chandler and Andersen (WCA) to partition 
the energy ( 5| into the short-ranged repulsive and the 
long-ranged attractive part [3], 
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The parameter A controls the strength of the attrac- 
tive solute-solvent interactions. 

We performed computer simulations using the 
Monte Carlo method at fixed pressure and tempera- 
ture for different values of the solute size R and inter- 
action strength A. The number N of solvent particles 
ranged from 864 to 8192 depending the size of the so- 
lute. The solute-solvent interaction ([6| was truncated 
and shifted at a distance r — R + 2.5(7. We calculate 
the average reduced solvent density at position r. 



g{r) = {p{r)) /p, 



(9) 



in an ensemble where the solute position is fixed at 
r — Q. Due to the isotropic nature of the solute this 
density depends only on the distance r = |r| from the 
solute center. 

We find that the reduced density g{r) has generally 
a well-defined first peak corresponding to the first sol- 
vation shell of the solute. The height of this peak. 
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FIG. 1: Magnitude of the first peak of tlie solvent density 
distribution function around liomogeneous Lennard-Jones 
solutes of radii Rja — 0.5, 0.7, 1, 2, 3, 5, plotted as a func- 
tion of tlie bare attraction strengtli A (top) and the min- 
imum of tlie solvent-solute potential (bottom). For large 
solute sizes and small values of A the reduced density dis- 
tribution function did not have a well-defined first peak, 
and no data is shown for those systems. 



denoted gmax, is shown in the top panel of Fig. [TJfor 
different values of R and A. This graph encompasses 
many of the important effects of our understanding of 
solvophobic and solvophilic solvation. 

First, for a fixed value of A less than 0.9, the solvent 
density near the solute decreases with increasing so- 
lute radius. It has been observed previously that the 
contact density near a hard sphere solute decreases 
from values above the bulk liquid density at small so- 
lute sizes to the equilibrium vapor density at large so- 
lute sizes, both for the Lennard-Jones solvent [10] and 
water [S] . The creation of a liquid-vapor interface near 
large enough solutes is the basis of our understanding 
of the hydrophobic effect pJLj . Our simulations show 
that this behavior persists for solutes with weak but 
finite attraction strengths A. 

Second, for any given solute size R, the density near 
the solute increases approximately linearly with the 
attraction strength A. This smooth response of the 
density to an attractive field over a wide range of in- 
teractions strengths shows that the solvent properties 
do not change significantly as A increases. In particu- 
lar it is consistent with Gaussian density fluctuations 
for all attraction strengths. 

Third, it is apparent that the slope 9(7max/9A, while 



FIG. 2: Space filling view of the solute clusters with parti- 
cle numbers n = 1, 13, 57, 135, 305. The clusters are spher- 
ical cuts of a hexagonally close-packed crystal. 



independend of A, depends on the radius R of the so- 
lute. In particular, there is a value A* « 0.9 of the in- 
teraction strength where the mean density itself seems 
to be independent of the solute size. A similar result is 
suggested by Ashbaugh and Paulaitis for the density 
near methane clusters in water |12j . This behavior can 
be easily understood within the framework introduced 
in Section |Tj The numerical value of A by itself is not 
an accurate parameter to describe the strength of the 
solute-solvent interactions. For the special case of the 
homogeneous Lennard-Jones solute, a much more re- 
liable indicator is the magnitude of the minimum in 
the potential energy, A|<i>niin(-R)|- The bottom panel 
of Fig. [l] shows that the density response is the same 
for all solute sizes if the control parameter is chosen 
to be the actual magnitude of the solute-solvent at- 
tractions. 

Our results for this simple model solute serve as a 
useful illustration for the theory outlined in Section [T] 
We now turn to more complicated solutes, where the 
lack of isotropy demands the analysis of more com- 
plicated density distribution functions. While these 
complications can obscur the underlying physics, we 
will see that the basic picture of solvation presented 
here remains correct. 



B. Cluster solutes in Lennard-Jones solvent 

We consider clusters of n particles solvated in the 
same Lennard-Jones fluid as before. The clusters 
are the model solutes of Ashbaugh and Paulaitis [T5] . 
Space flUing views are shown in Figure [2j The clusters 
are composed of n — 1,13,57,135, and 305 hexago- 
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nally close packed particles with nearest neighbor dis- 
tance d = 2^/^(T. The solute interacts with the solvent 
molecules with the potential energy 
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where is the position of the ath cluster particle, 
Ri is the position of the ith solvent particle, and 



wx{r) = ^0(7") + Xwi{r). 



(11) 



Here, WQ{r) and wi{r) are, respectively, the repul- 
sive branch and the attractive branch of the Lennard- 
Jones potential 3., 








if r < 2i/Vs 
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if r < 21/6(7, 
if r > 2i/6cr, 



(13) 

We choose as — a and eg = e so that the interaction 
between a cluster particle and a solvent particle is the 
same as that between two solvent particles at A = 1. 

The n — 1 solute is spherically symmetric, and 
the mean solvent density near it is described by the 
isotropic radial distribution function g(r) as before. 
For the clusters with n > 1, however, specification of 
the mean solvent density requires many more variables 
than r, and it is natural to consider projections of the 
mean density. One such projection is the so-called 
proximal distribution gproxir), defined as the reduced 
density averaged over the surface [T^J [13] 
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where 



(r) = min{|r - r„| : a = 1, . . . ,n} (15) 



is the distance from a point r to the nearest clus- 
ter particle. While explicitly a function of only one 
length, r, it is implicitly a multi-point distribution 
function [M], except in the case n — 1. For that spe- 
cial case (7prox('') is the standard radial distribution 
function g{r). 

Figure [3] illustrates the behavior of the proximal 
distribution function for the model solutes in the 
Lennard- Jones solvent. Ashbaugh and Paulaitis pro- 
posed the height 5,nax of the first peak in gproxif) 
as a meaningful measure of the solvent density near 
arbitrary solutes. In the top panel we show that 
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FIG. 3: Proximal solvent density distribution function 
around clusters of Lennard-Jones particles immersed in 
the Lennard-Jones fluid, (top) Height of the flrst peak 
in gprox(r) for different cluster sizes n and attraction 
strengths A. (bottom) Proximal distribution functions for 
A = A* = 0.5. 



the behavior of this observable is qualitatively the 
same as that shown in Fig. [l] for isotropic so- 
lutes. The response of the mean solvent density to 
changes in solute-solvent attraction is proportional to 
that change. This result is as expected according to 
Eq. ([2]) . Another important point is that this propor- 
tionality leads to a coincidence of curves at A = A*, 
where in this case A* « 0.5. The existence of a coin- 
cident value of A is anticipated in Eq. 

Ashbaugh and Paulaitis found in their study of 
methane clusters solvated in water that at A* not 
only the heights of the first peak in 5prox('") coincide 
for different cluster sizes, but that the proximal dis- 
tribution functions themselves are very similar. The 
bottom panel of Fig. |3] establishes a similar result for 
our model system. From this behavior, Ashbaugh and 
Paulaitis conclude that the solvation environment of 
a monomeric methane particle is essentially the same 
as that of methane confined to a cluster of arbitrary 
size. Fig. |4]shows that this conclusion is not justified. 
The contour lines of constant density do not coincide 
with the surfaces of constant proximal distance for the 
n = 135 cluster. In particular, the crevices of the so- 
lute, where the attraction strength is large, contain 
a much higher solvent density than the proximal dis- 
tribution function suggests. No such behavior exists 
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FIG. 4: Reduced solvent density around the n = 135 
solute cluster immersed in the Lennard-Jones fluid at A = 
A* = 0.5. Shown is a slice through the solute center of 
width O.lo" with resolution 0.05a, demonstrating that the 
solvent density does not follow the shape of the solute 
surface. 



for the monomeric solute, where the density profile is 
isotropic. 

As the strength of the attactive solute-solvent in- 
teraction increases, the clusters lose their hydrophobic 
character until eventually they become hydrophilic. 
To quantify this process we calculate the solvation free 
energy, or excess chemical potential, as a function of 
A, which is given by 



A/ZA =Afio+ dA' {Wi - Wo); 



where 



(16) 



indicates ensemble average with solute- 



solvent attractive force strength A. 

Solvation free energies due to excluded volume are 
approximately proportional to the excluded volume 
for small solutes. Further, at large excluded volumes 
with small surface curvature and no attractive solvent- 
solute forces, solvent density is depleted near the so- 
lute. Therefore, for the solutes we have considered, 
A/xq can be well approximated by the solvation free 
energy for a hard sphere with the same excluded vol- 
ume as the solute clusters [TU]. Here we define the 
excluded volume as those spatial regions where the 
solute-induced field is greater than the thermal en- 
ergy /cbT. 

The result of this calculation is shown in Fig. [5j 
Over the range of interaction strengths considered, 
the solvation free energies change from strongly pos- 
itive to negative. In other words, at large values of 
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FIG. 5: Excess solvation energy of solute clusters in the 
Lennard-Jones fluid. The symbols are calculated from 
simulation using Eq. (161, and the lines are quadratic 
fits through the data. 



A the clusters are not solvophobic. This transition 
occurs smoothly and regularly. The shape of these 
curves is understood, because the attractive energy 
is approximately proportional to A times an integral 
over {p{r))^, and the latter is a linear function of A, 
as noted in Eq. 



C. Cluster solutes in SPC water 

In another series of calculations, the solvent is the 
SPC model for liquid water [15J, and the solutes are 
clusters of methane molecules in the united atom ap- 
proximation of Ref. 16. As in Ref. [T^ these clusters 
are hexagonally close-packed with a particle spacing 
of c? = 4.I9A. The methane molecules interact with 
the water oxygens through a Lennard-Jones poten- 
tial with cr = 3.448A and e = 0.895kJ/mol, which 
was truncated and shifted at a distance of 9A. As 
before we use the WCA partition of the interaction 
potential to obtain control over the strength of the 
attractions, given by the parameter A. Electrostatic 
interactions were calculated using Ewald summation. 
Molecular dynamics simulations were performed us- 
ing the LAMMPS program [171 CH] at temperature 
T = 298K and a box size chosen to give the correct 
bulk density [H]. 

As noted above, Ashbaugh and Paulaitis [T2] found 
that the proximal solvent density distribution func- 
tions are very similar for different cluster sizes at the 
specific attraction strength A* w 1, from which they 
concluded that the solvent density near a cluster of 
solute particles is basically the same as that near a 
single solute particle. In Fig. |6]we show the reduced 
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FIG. 6: Reduced water density around the n — 135 solute 
cluster at A = A* = 1. Shown is a slice through the solute 
center of width 0.5A with resolution 0.25A. To a good 
approximation the solvent density follows the shape of the 
solute cluster. 
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FIG. 7: Excess solvation energy of methane clusters in 
water. The sy mbo ls are computed from simulation with 
the aid of Eq. ( 16 1 and known solvation free energies of 



hard sphere solutes in water [27]. The lines are quadratic 
fits to the data. 



III. DISCUSSION 



water density near the n = 135 solute with an in- 
teraction strength A = A*. For this specific system 
the solvent density is indeed approximately constant 
along contours of constant proximal distance, so that 



(p(r)) w pgprox(Vox('')) 



(17) 



is a rather good approximation for the range of pa- 
rameters considered. However, this agreement is not 
exact, as deviations are visible in the crevices. In these 
regions adjacent to the cluster, the solvent interacts 
with many solute particles, and the attractive poten- 
tial well reaches values as deep as — 5kJ/mol. This 
value is to be compared with the average energy of 
a water molecule in the bulk fluid, which is approxi- 
mately 40kJ/mol 112]. Thus, the cluster crevices are 
approaching hydrophilic character and can compen- 
sate for a fraction of the average water-water interac- 
tion energy. 

Fig. |7] shows the solvation free energy of the 
methane clusters in water, calculated using ( 16 1 as 



before. Again the clusters show a smooth transition 
from hydrophobic solvation at small values of A to hy- 
drophilic solvation at large attraction strengths. For 
this particular system, the crossover between these 
two regimes occurs at values of A greater than A*. 



The results of the previous section show that the 
role of solvent-solute attractions on solvation is sig- 
nificant but simple. In particular, while fluctuations 
of solvent density are affected by the shape and size 
of the solute and the thermodynamic state of the sol- 
vent, the mean solvent density changes only linearly 
with changes in the strength of adhesion between so- 
lute and solvent. If, for example, a solute is small, the 
hquid solvent will be constrained to be outside the 
solute's excluded volume but otherwise behave as if 
it were bulk. This conception is precisely the physics 
that underlies the Pratt-Chandler theory of hydropho- 
bicity [7| 120] and its contemporary variant [5] . On the 
other hand, if the solute's girth is such to exclude a 
volume with large low-curvature surface area, fluctu- 
ations of the solvent are changed to reflect the forma- 
tion of an interface jlj [TTJ [21] ■ To the extent there is 
no solvent-solute adhesion, the interface is essentially 
the same as that between liquid and vapor |22j . The 
role of attractions between the solute and solvent is 
not to change the fluctuations, but to affect the aver- 
age solvent density and the position of this interface. 
For strong enough attractions, the effect is to pin the 
interface at speciflc positions, as seen in Fig. [4] and 
to a lesser extent in Fig. |6] 

This picture of the role of attractions in the con- 
text of hydrophobic solvation has been discussed be- 
fore [Hi E]. Simulation results of Werder et al. on 
the hydration of several different models of graphite 
in water p3i can be viewed as illustrations of that role. 
Namely, Werder et al. establish linear trends for hy- 
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dration energies as functions of solvent-solute attrac- 
tive energy strength. The relative range of parameters 
is smaller in that work than in ours, probably explain- 
ing why curvature in the linear trends is not apparent. 
Along with demonstrating regularity and simplicity, 
Werder et al. also show that the trends are significant 
in that modest potential parameter changes result in 
models of "graphite" that change from hydrophobic 
to hydrophilic. In one model, Werder et al. consider 
a Lennard-Jones potential between graphite carbon 
atoms and water oxygen atoms with energy and length 
parameters eco = 0-39 kJ/mol and ctco = 0-32 nm. 
They establish that in this model, line tension is pos- 
itive and the contact angle between water and a pla- 
nar graphite surface is finite. In contrast, a similar 
model of a presumed hydrophobic surface, but with 
£co — 0-48 kJ/mol and aco = 0-33 nm yields a 
negative surface tension between water and the sur- 
face [23] |5S] . This latter surface is therefore not hy- 
drophobic. 

An important consequence of the simplicity we em- 



phasize regards the entropy of solvation. To the extent 
that density fluctuations obey Gaussian statistics (i.e., 
the mean density responds linearly to forces of adhe- 
sion) , the entropy of solvation is independent of adhe- 
sive energy strength. This idea is used in Ref.[2B]to re- 
place the hydration entropy of a protein with that of a 
hard sphere occupying the same volume, even though 
the surface of the folded protein is largely hydrophilic. 
By demonstrating linear trends for the mean density, 
our work here provides explicit justification for this 
replacement. 
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In Ref. 1231 Werder et al. do not report the oil-water 
surface tension for these particular parameters, but 
Choudhury and Pettitt, in effect, do so in Ref. 1241 
Specifically, for carbon-carbon potential parameters, 
the authors of [53] cite ecc = 0.36 kJ/mol and acc = 
0.34 nm, and combine these with oxygen-oxygen pa- 
rameters eoo ~ 0.65 kJ/mol and aoo = 0.32 nm, 
using geometrical and arithmetic means, respectively, 
to obtain their carbon-oxygen parameters. They com- 



pute a potential of mean force by molecular simulation 
for a pair of plates formed from these "carbon" parti- 
cles. Figure, la of their paper shows that the solvent 
contribution to the free energy for dissociated plates 
is lower than that of associated plates. The plates 
are reasonably large, so that the solvation energy is 
predominantly proportional to the solute surface area 
exposed to the solvent. Choudhury and Pcttitt's re- 
sult for the solvent contribution to the potential of 
mean force therefore shows that the interactions used 
in their study produce a negative solvation surface 
free energy for the hydrated plate. 
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